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Abstract 

Powder-snow avalanches are violent natural disasters which represent a major risk 
for infrastructures and populations in mountain regions. In this study we present 
a novel model for the simulation of avalanches in the aerosol regime. The second 
scope of this study is to get more insight into the interaction process between an 
avalanche and a rigid obstacle. An incompressible model of two miscible fluids can be 
successfully employed in this type of problems. We allow for mass diffusion between 
two phases according to the Fick's law. The governing equations are discretized 
with a contemporary fully implicit finite volume scheme. The solver is able to deal 
with arbitrary density ratios. Encouraging numerical results are presented. Volume 
fraction, velocity and pressure fields are presented and discussed. Finally we point 
out how this methodology can be used for practical problems. 
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Figure 1. Two illustrations of powder-snow avalanche flows. 
1 Introduction 

Snow avalanches are commonly defined as abrupt and rapid gravity-driven 
flows of snow down a mountainside, often mixed with air, water and some- 
times debris (see Figure 1). Avalanches are physical phenomena of great in- 
terest, mainly because they represent a big risk for those who live or visit 
areas where this natural disaster can occur. During last several decades, the 
risk has increased due to important recreational and construction activities in 
high altitude areas. We remember recent events at Val dTsere in 1970 and in 
Nothern Alps in 1999 [AncOl]. 

The avalanches arise from an instability in a pile of granular material like 
sand or snow [TM08]. The destabilization phase of an avalanche life is still a 
challenging problem. There are many factors which influence the release pro- 
cess. One can recall snowpack structure, liquid contain, shape and curvature 
of starting zone and many others [AncOl]. In this study we focus especially on 
sliding and stopping phases. 

The serious research work on this natural phenomenon was preceded by the 
creation of scientific nivology at the end of the XlXth century. Among the 
pioneers we can mention Johann Coaz (swiss engineer) [Coa81] and P. Mougin 
(french forest engineer, author of the first avalanche model using an analogy 
with a sliding block) [Mou22,Mou31]. The work of P. Mougin was ignored until 
the 1950s when A. Voellmy developed a similar model [Voe55]. Its somehow 
improved versions are still used by engineers nowadays. 

We would like to point out here important contributions of the Soviet school by 
S.S. Grigorian, M.E. Eglit, A.G. Kulikovskiy, Y.L. Yakimov and many others 



[GEY67,KE73,BE73,BKK + 75,KS77b,Egl83,Egl91,BL98,Egl98,BEN02]. They 
were at the origin of all modern avalanche models used nowadays in the en- 
gineering practice and, sometimes, in scientific research. Their works were 
mainly devoted to the derivation and comprehension of mathematical models 
while occidental scientists essentially looked for quantitative results. 

Conventionally we can divide all avalanches in two idealized types of motion: 
flowing and powder-snow avalanches. A flowing avalanche is characterized by 
a high-density core ranging from 100 to 500 kg/m 3 and consists of various 
types of snow: pasty, granular, slush, etc. The flow depth is typically about a 
few meters which is much smaller than the horizontal extent. This argument is 
often used to justify numerous depth-integrated models of the Savage-Hutter 
type 1 [SH89,SH91]. These avalanches can cause extensive damage because of 
the important snow masses involved in the flow in spite of their low speed. 

On the other hand, powder-snow avalanches are large-scale turbidity currents 
descending slopes at high velocities [RH04]. They seriously differ from flowing 
avalanches. These clouds can reach 100 m in height and very high front veloc- 
ities of the order of 100 m/s. They grow continously and the average density 
is fairly low (from 4 to 25 kg/m 3 ). These spectacular avalanches (see Figure 

1) occur only under certain conditions (after abundant fresh snowfalls, cold, 
dry and weakly cohesive snow on strong slopes) and they produce a devastat- 
ing pressure wave which breaks the trees, buildings, tears off the roofs, etc. 
During the propagation stage, they are able to cross the valleys and even to 
climb up on the opposite slope. Hence, measurements by intrusive probes are 
almost impossible. Avalanches in aerosol are not very frequent events in Alps 
but in the same time we cannot say they are very seldom. In the technical 
literature there is an opinion that an avalanche in aerosol is less destructive 
than a flowing one since the transported mass is much smaller. Nevertheless, 
recent events of the winter 1999 in Switzerland, Austria and France revealed 
the important destructive potential of the powder-snow avalanches (see Figure 

2) . 

Recently several systematic measurements campaigns in situ were conducted 
in Norway, Switzerland and Japan [MS84,NTK90,NMKI93,NSKL95,DGBA00]. 
Researchers shed some light on the internal structure of big avalanches. More 
precisely, they show that there exists a dense part of the avalanche which re- 
mains permanently in contact with the bed. This dense core is covered by the 
aerosol suspension of snow particles in the air. From these results it follows 
that mentioned above two types of avalanches may coexist in nature and pro- 
posed above classification is rather conventional. Perhaps, future studies will 
perform a coupling between the dense core and powder-snow envelope in the 



In hydrodynamics and hydraulics this type of modeling is also known as shallow 
water or Saint- Venant equations [dSV71]. 
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Figure 2. Protecting wall in armed concrete at Taconnaz (Haute-Savoie, France) 
destroyed by the powder-snow avalanche of the 11 February 1999. The height is 7 
m and the thickness is 1.5 m (Photo by C. Ancey). 

spirit of SAMOS code [SZ04]. 

Let us review various existing approaches to the mathematical modeling of 
snow avalanches. Generally, we have two big classes of mathematical mod- 
els: probabilistic and deterministic. In the present article we deal with a 
deterministic model and we refer to the works of K. Lied and D.M. Mc- 
Clung [LB80,ML87,McC00,McC01] for more information on statistical ap- 
proaches to avalanche modeling. Deterministic models can be further divided 
into continuous and discrete ones depending whether the material under con- 
sideration can be approximated as a continuum medium or not. For the 
review and some recent results on dense granular flows we refer to follow- 
ing works [NKMN98,MN01,Raj02b,Raj02a,Raj05] and the references therein. 
Some promising results were obtained with discrete models based on cellular 
automata [DGRS + 99,ADGM+00,DSI06]. 

The first contemporary avalanche models appeared in 1970 by soviet scientists 
Kulikovskiy and Sveshnikova [KS77b,BL98]. Later, their idea was exploited 
by Beghin [Beg79,BB83,B091] and others [HTD77,FP90,AU99,Anc04,RH04]. 
We call this type of modeling OD-models since the avalanche is assimilated 
to semi-elliptic cloud with variable in time volume V(t), momentum (pU)(t), 
etc. All quantities of interest are assimilated to the center of mass and their 
dynamics is governed by conservation laws expressed as Ordinary Differential 
Equations (ODE). 

On the next complexity level we have various depth-integrated models. The 
governing equations are of Shallow Water (or Saint- Venant [dSV71]) type. 
In general, they are derived by depth averaging process or some asymptotic 
expansion procedure from complete set of equations. Thus, a physical 3D (or 
2D) problem results in a 2D model (ID correspondingly). From computational 
point of view these models are very affordable even for desktop computers. On 
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the other hand, they provide us very approximative flow structure, especially 
in the vertical direction. 

In France, G. Brugnot and R. Pochat were among the pioneers [BP81] while 
in Soviet Union this direction was explored in the beginning of 1970 by N.S. 
Bakhvalov et M.E. Eglit [BE73,BKK + 75]. Currently, each country concerned 
with the avalanche hazard, has its own code based on this type of equations. 
Probably the most representative model is that developed by S. Savage and 
K. Hutter [SH89]. Nowadays this set of equations is generally referred as the 
Savage-Hutter model. We refer to [Hop83,Hut96] as general good reviews of 
existing theoretical models and laboratory experiments. 

This approach was further developed by incorporating more complex rheolo- 
gies and friction laws [GWH98,Hut91,HG93,MCVB+03,HWP05,FNBB + 08]. 
To conclude on this part of our review, we have to say that this modeling is 
more relevant to the flowing avalanche regime which is characterized by the 
small ratio of the flow depth h to the horizontal extent £ (i.e. k <C 1). 

Finally, we come to the so-called two-fluid (or two-phase) models. In this 
paradigm both phases are resolved and, a priori, no assumption is made on 
the shallowness of the flow under consideration. Another advantage consists in 
fact that efforts exerted by the ambient air on the sliding mass are naturally 
taken into account. From computational point of view, these models are the 
most expensive [NG98,ESH04,Eti04,EHS05]. In the same time, they offer quite 
complete information on the flow structure. 

As it follows from the title, in this study we are mainly concerned with powder- 
snow avalanches. We would like to underline that our modeling paradigm al- 
lows for taking into account the dense core at the first approximation order. 
The density is completely determined by the snow volume fraction distribu- 
tion. This parameter can be used to introduce a stratification in the initial 
condition, for example. Otherwise, it will happen automatically due to mixing 
and sedimentation processes in the flow, provided we follow its evolution for 
sufficiently long time. 

We retained a simple incompressible two-fluid model which is described in 
detail in Section 2. Two phases are allowed to interpenetrate, forming a mixing 
zone in the vicinity of the interface. We make no Boussinesq-type hypothesis 
[Anc04] about the density ratio. Moreover, our solver is robust and is able 
to deal with high density ratios (we tested up to 3000). In nature, the flow 
under consideration is obviously highly turbulent but at the present stage we 
do not incorporate any explicit turbulence modeling beyond resolved scales. 
In the present study we focus mainly on the physically sound description of 
incompressible highly inhomogeneous two-phase flows. 

Good understanding of these natural phenomena can improve the risk assess- 
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ment of such natural hazards. Large scale experiments are not feasible 2 . Field 
measurements during the event are very hazardous and the events are rare. 
Laboratory experiments on powder-snow avalanches are essentially limited to 
Boussinesq clouds 3 [Kel95,NBNBH02,Pri03,NB03,PNBNF04] while it is not 
the case in nature. Thus, we do not really know how these results apply to real- 
world events. Fortunately, some progress has been made recently to remedy 
this situation [TM08]. 

Numerical simulations of avalanches provide useful information on the dynam- 
ics of these flows. Computer experiments may become in term the main tool 
in testing various situations. Experiments in silico should be complementary 
to those in situ or in laboratory. Direct Numerical Simulations (DNS) provide 
complete information about all flow quantities of interest such as the local den- 
sity, the velocity field variations, the dynamic pressure and the energy. Recall 
that this information is not easily accessible by means of measurements. 

The present paper is organized as follows. In Section 2 we present the govern- 
ing equations and some constitutive relations. Special attention is payed to the 
presence of strong density gradients in the flow and to the kinetic energy bal- 
ance of the resulting system. The next Section 3 contains a brief description of 
the numerical methods and numerous computation results are presented. Fi- 
nally, this paper is ended by outlining main conclusions and a few perspectives 
for future studies (see Section 4). 



2 Mathematical model 

In the present study we assume that an avalanche is a two-fluid flow formed by 
air and snow particles in suspension. The whole system moves under the force 
of gravity. For simplicity we assume that the mixture is a Newtonian fluid. 
The last assumption is not so restrictive as it can appear. The flow under 
consideration is such that the Reynolds number is very high (Re ~ 10 9 ). 
Therefore, the transient behaviour is essentially governed by the convective 
terms and not by the fluid rheology. On the contrary, the rheology is very 
important in the flowing regime. 



2 However, we would like to notice that there are two experimental sites: one in 
Switzerland (Sion Valley) and another one in France (Col d'Ornon). Unfortunately, 
field measurements provide very limited information on the flow structure [DGA01], 
even if some a posteriori analysis may be beneficial [TM07]. 

3 The Boussinesq regime corresponds to the situation when p+ ~+ <C 1 where 
are densities of the heavy and light fluids respectively. This asymptotics allows to 
introduce the so-called Boussinesq approximation. 



6 




Figure 3. An elementary fluid volume dVl occupied by two phases. 

In two-fluid flows it is natural to operate with the so-called volume fractions 
[Ish75,TK96,TKP99,GKC01]. Consider an elementary fluid volume dVt sur- 
rounding an interior point P G dVt. Let us assume that the first fluid occupies 
volume dVt\ C dQ and the second the volume dQ 2 ^ dfl (see Figure 3) such 
that 

\<Kl\ = |dfii| + \dtt 2 \. (1) 
The volume fraction of the fluid i — 1, 2 in the point P is defined as 



UP) := lim 



idni->o \dQ\ 

PGdQ 

From relation (1) it is obvious that <fri{P) + 02 (-P) = 1, for any point P in 
the fluid domain. Henceforth, it is sufficient to retain only the heavy fluid 
volume fraction 0i, for example, which will be denoted just by 0, for the sake 
of simplicity. 

If we assume that constant densities and kinematic viscosities of the heavy 
and light fluids are respectively p ± and u ± , the mixture density p and the 
dynamic viscosity p are defined as follows: 

p = # + + (l-0)p-, p = <Pp + v + + (l-(j>)p-v-. (2) 

In the following we assume that p + ^ p~ , otherwise the two-fluid modeling 
does not make much sense. 

After some simple algebraic computations, the mixture dynamic viscosity p 
can be expressed in terms of the mixture density p as: 

p = p + up, (3) 

where po has the dimension of the dynamic viscosity [kg/(m ■ s)] and v scales 
with the kinematic one [m 2 / s\. These coefficients are related to the densities 
and kinematic viscosities of constitutive fluids in this way: 

v~p~ p + — v + p + p~ _ u + p + — v~p~ 



Po := t , v :-- 



Compressible Navier-Stokes equations with viscosities of the form (3) were 
studied mathematically in [Sy05]. 
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Remark 1 In view of our applications, the snow kinematic viscosity u + can 
be parameterized as a function of temperature T and snow density p + according 
to [DAH+07J: 

P + 

where p = 3.6 x 10 6 TV • s ■ m~ 2 , a = 0.08 K~ l , (5 = 0.021 m 3 /kg. 

Now we can state the governing equations of our physical problem. In this 
study we assume that both phases are constrained to have the same velocity 
variable u. This assumption is not very restrictive. The formal justification 
of single- velocity two-phase models can be found in [MG05,MDG09,MDG10]. 
Also, this type models have already been successfully applied to a number of 
practical problems [Dut07,DDG08b,BPB09,DDG10,DDG08a]. 

The flow is assumed to be isentropic. The mass and momentum conservation 
equations have the classical form: 

dtp + V • (pu) = 0, (4) 

d t (pu) + V ■ (pu <g> u) + Vp = V ■ (2pB(u)) + pg, (5) 

where g is the acceleration due to gravity, p is the hydrodynamic pressure and 
D(m) = |(Vu + (Vw)'j is the strain rate tensor. 

The fluid mixing is taken into account by Fick's type law [Fic55a,Fic55b] 
resulting in the following quasi-compressible equation: 

V-u = -V- (/cVlogp), (6) 

where the coefficient k has the dimension of kinematic viscosity and will be 
defined below. For the moment, we assume k to be constant, consequently, the 
right-hand side of equation (6) can be equivalently rewritten as — kA log p. 

Remark 2 Examples of closures similar to (6) of the form 

V • u = ±V • (V0(p)) 

may be found m [Gra55,KS77a,BadV83,Maj84,AKM90,FS01] (the sign ± de- 
pends on monotonicity properties of the function <j)(p) ). It may lead to models 

of low Mach number combustion [Maj84] with cj)(p) = —, salt or pollutant mo- 
tion in a shallow layer [FS01] and other practical situations involving highly 
mhomogeneous fluids [BadVSV82,BadV83,BES07j. 

Remark 3 Since we model an avalanche propagation along a sloping solid 
boundary, we take the gravity acceleration in the following form: 

g = (#sin0, -gcos9), 
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where 9 is the slope of the hill and usually g := \g\ = 9.81 m/s 2 . 
2. 1 Model based on fluid volume velocity 

In this section we transform governing equations (4), (5) and (6) to oper- 
ate with more representative physical variables giving also a more convenient 
mathematical form. 

Namely, we introduce the new velocity variable defined as 

v := u + KVlogp, (7) 

which is sometimes referred in the literature as the mean volume velocity (cf. 
[FS01]) or the fluid volume velocity (cf. [Bre05a]). In this study we will retain 
the last term. It has been pointed out recently [Bre05b,Bre06] that the fluid 
volume velocity v is more pertinent for flows involving high density gradients. 
Now, let us rewrite the system (4), (5) and (6) in terms of new variable v. 

First of all, from Fick's law (6) it follows immediately that the flow is incom- 
pressible within the fluid volume velocity v: 

V • v = 0. 

The mass conservation equation (4) takes the following simple form: 

dtp + V ■ (pv) = V ■ (k P V logp). 

By taking into account the fact that k is constant and the field v is divergence- 
free (V • (pv) = v ■ Vp), we can rewrite the last equation as: 

d t p + v ■ Vp = i%Ap. (8) 

The latter equation is of parabolic type. The diffusive term comes from the 
Fick's law governing the mixing process between two fluids. In the case when 
k = 0, the initial sharp interface between two phases would be simply advected 
by the velocity field v (coinciding with u, when k = 0), thus, preventing any 
mixing. Henceforth, we consider the case k > 0. 

Remark 4 The mass conservation equation (8) can be equivalently rewritten 
in terms of the volume fraction <fi using the mixture density representation (2): 

d t (j) + v ■ V0 = kA0. (9) 

Finally, we have to transform the momentum conservation equation (5) ac- 
cording to the change of the velocity variable (7). This operation will require 
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several computations briefly presented below. For convenience, we work with 
equation (5) recasted in non-conservative form using the mass conservation 
(4): 

pd t u + p(u ■ V)w + Vp = pg + V • (2/iD(m)). (10) 

The three terms pd t u, p(u ■ V)u and V • (2pD(w)) involving u have to be 
rewritten: 



pd t u = pd t v - Kpd t (—\ (11) 
v P 1 

p(u ■ V)u = p(v ■ V)v — Kp(V log p ■ W)v — Kp(v ■ V)V log p (12) 

+ K 2 p(Vlogp- V)Vlogp, (13) 

V • (2pB(u)) = V • (2fjB(v)) - kV • (2pVVlogp). (14) 

Vp 

In order to obtain the evolution equation for the quantity arising in (11), 

P 

we use the mass conservation equation (8): 

pd t (—)+ V(V- Vp) - (v-Vp)— = KAVp - kA P ^. 

v p ' v ' p p 

Consequently, relation (11) can be rewritten using the last result: 

pd t u = pd t v + kV (v ■ V p) - n(v-Vp)— + k 2 — Ap - n 2 AVp. 

P P 



After all these developments, the momentum conservation equation (10) be- 
comes: 

pd t v + p(v- V)t7+ Vp + kV(v- Vp) - k(v- Vp)— + k 2 — Ap-n 2 AVp 

V v ^ V ^ ' 

(I) (II) 

— Kp(V log p • V)v — Kp(v ■ V) V log p + K 2 p(V log p • V) V log p = 

pg + V • (2pD({/)) - kV • (2pVVlogp). 

One can remark that terms in groups (I) and (II) can be simplified to give 
K^VvVp and k 2 V ■ \ <g> pj correspondingly: 

pd t v + p(v ■ V)v + Vp + n l VvVp - nVpVv - n 2 AVp + n 2 V ■ f — <g> p 

v P 



(*) 

pg + V- (2pB{v)) -2vkV ■ (pVVlogp) -2/tp V • (VVlogp) . (15) 

v v ' v „ ' 

(*) (**) 
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In order to obtain last two terms on the right-hand side of (15), expression 
(3) for the mixture dynamic viscosity p was used. 

However, equation (15) can be further simplified if we make a special choice 
for the constant k arising in the Fick's law (6). More specifically, we take 
k = 2i>, where v is the mixture kinematic viscosity defined in equation (3). 
For this choice of Fick's diffusion coefficient k, the terms marked with (*) in 
(15) disappear, since it is straightforward to check the following differential 
identity: 

V ■ (pVVlogp) = AVp - V ■ (— ® p). 

v P ' 

Finally, the term (**) in equation (15) can be written as a gradient: 

V ■ (VVlogp) = VAlogp. 
Consequently, it can be incorporated in the definition of the pressure: 

7r(x, t) := p(x, t) + 4vp A log p. (16) 

If we summarize developments made above, we can state the governing equa- 
tions of the proposed model: 

V-t7=0, (17) 
d t p + v- Vp = 2uAp, (18) 
pd t v + p(v- V)v + Vtt + 2v t VvVp - 2u\7pS7v = pg + V ■ (2/iD(t7)) (19) 

Equations (17), (18) and (19) have to be completed by appropriate initial and 
boundary conditions to form a well-posed problem. 

Initially, the velocity, pressure and density (or equivalently, volume fraction) 
fields have to be imposed. Concerning boundary conditions, the usual no-slip 
condition v — can be imposed. Another possibility is to impose the following 
partial slip condition: 

v-n = 0, ((1 - a)v + a(B(v) • n)) •? = 0, (20) 

where f is the tangent vector to the boundary and a G [0, 1] is the friction 
parameter. We have to say that the latter is known to be physically more 
relevant [Eti04]. 

If it is required by a computational algorithm as in our case, for example, 
governing equations can be recast in the conservative form. Using the incom- 
pressibility condition (17), the mass conservation equation (18) becomes: 

d t p + V- (pv) = V- (2PVp). 
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If we sum up equation (19) with mass conservation (18) multiplied by v, we will 
obtain the following equation with advective terms written in the conservative 
form: 

d t (pv) + V • (pv <g> v) + Vvr + 2v t Vv Vp - 2z/VpW = 

pg + 2vvAp + V ■ (2pB(v)). (21) 

The dimensionless form and discretization of these equations will be discussed 
briefly below in Sections 2.3 and 3 correspondingly 



2.2 Kinetic energy evolution 



In this section we would like to derive an integral identity which describes the 
kinetic energy evolution associated to the system (17) - (19). Throughout all 
developments in this section we assume the no-slip condition v = on the 
fluid domain Q boundary dfl. The same result will hold if we assume periodic 
boundaries or appropriate decay conditions at the infinity. 

\v\ 2 

First, we multiply the mass conservation equation (18) by , momentum 

conservation (19) by v, sum them up and integrate over the fluid domain Q. 
After a few integrations by part, using the incompressibility (17) and boundary 
conditions, we obtain the following identity: 

d t J P^- dx = J pg-vdx + J 2vAp^—^- dx + J 2u(VpVv)vdx + 
n n n n 

" v ' 

CO 

V • (2upB{v))vdx - J 2v( t VvVp)vdx + J V ■ {2p ^{v))vdx 



(ii) (in) 
Two terms from the group (I) cancel each other after one integration by parts: 

/I I C I ^/ 1 v 

2uAp l -^-dx = - J 2uVp- v(^-)df = -J 2u(VpVv)vdx. 

a n n 

The group (II) terms can be also transformed to give: 

-J 2up\3(v)\ 2 dx + J 2upWv: Vvdx = - J 2up\A{v)\ 2 dx, 



where A(v) := |fVu — *Vuj is the antisymmetric part of the velocity gradient 
Vv and A : B = a^bij is the usual contracted product of two square matrices 
A = (aij)i<i,j<n, B = (bij)i<ij< n . 
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Finally, the term (III) similarly can be transformed as: 

J V ■ (2p 3(v))vdx = - J 2p \B(v)\ 2 dx 
n n 

Thus, we come to the following integral identity governing the kinetic energy 
evolution associated to the system (17) - (19) (in closed conditions): 

9t I P 2 dS= f Pd'Vdx- ffyoP{v)\ 2 dx- J 2v p\k{v)\ 2 dx . (22) 
n n n n 

v ' S v ' 

(a) (b) 

In this integral equality each term has a specific physical sense. The term on 
the left-hand side represents the rate of kinetic energy change per a unit of 
time. On the right-hand side the first term is the work done by the gravity 
force, thus, describing the energy input into the system. Two last terms repre- 
sent various dissipative processes. The first term (a) can be seemingly ascribed 
to the dissipation due to viscous forces, while the second one (b) comes from 
mixing processes between two phases. 

Despite the fact that governing equations (17) - (19) are more complicated 
than classical Navier-Stokes equations, due to the judicious choice of Fick's 
constant k we were able to get a kinetic energy balance in simple and physically 
sound form (22). 

2.3 Dimensional analysis 

In this section we perform a dimensional analysis of the governing equations 
(17) - (19) in order to reveal important scaling parameters. Henceforth, starred 
variables denote dimensional quantities throughout this section. 

The initial avalanche height ho and the heavy fluid density p + are chosen to be 
the characteristic length and density scales correspondingly. The velocity field 
is adimensionalized by a typical flow speed Uq. Finally, from characteristic 
length and velocity, it is straightforward to deduce the time scale. Conse- 
quently, the scaling for the independent variables is 

-mi i -> ,* ho , 
x = hox, t = — t, 
u 

and dimensionless dependent variables p, u, p are introduced in this way: 

p* = p + p, U* = Uq U, 7T* = p + Ug 7T, /!* = p 4 D pL 
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The governing equations (17) - (19) in dimensionless form become: 



V-iT=0, 

d t p + v ■ Vp = — Ap, 
Re 

pd t v + p(v-V)v + Vit + -|-*WVp - -J-VpW= ip# + ^-V • (pO(*7)) 

Re Re rr Re 

This procedure reveals two important scaling parameters - the Reynolds num- 
ber Re [Rey83] and the Froude number Fr [B28] which are defined as 

Re:=^, Fr, U ° 



2u VPo 

Recall that the Reynolds number Re gives a measure of the ratio of inertial 
forces to viscous forces. In practice this number characterizes different flow 
regimes, such as laminar or turbulent flow. The Froude number is a ratio of 
inertia and gravitational forces. It is a hydrodynamic equivalent of the Mach 
number. For free-surface flows it specifies the nature of the flow (subcritical 
or supercritical) [DVB02]. 

Remark 5 We would like to discuss the difference between the physical pres- 
sure field p(x,t) and the modified pressure 7r(x,t). If we turn to dimensionless 
variables in equation (16) , we will obtain the following relation: 

7r(£ , t) = p(x, t) + —2 A log p. 

Re 

Taking into account the typical values of the Reynolds number in our appli- 
cations, we can conclude that this difference is completely negligible, at least 
from the practical point of view. 

In fact, there are two additional scaling parameters S and A hidden in defini- 
tions of the density p and dynamic viscosity of the mixture: 

p:=4 = < ^+( 1 -^ 5:= ^T> 
p+ p+ 

p + V 1 — OA V+ 

where we substituted the following representation of the Fick's coefficient: 

,1-5X 

v = ./+— . (23) 



Actually, the densities ratio S can be related to the well-known Atwood number 
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At (cf. [GGL+01,LJ05]): 



At := ^ = \± . 

P + + p~ 1 + 6 

The powder-snow avalanche regime is characterized by very high values of 
Reynolds number Re and low density ratios 5 [MS93,Anc03]: 

Re ~ 10 6 , 0.05 < 5 < 0.25. 



Thus, the flow is clearly turbulent. Nevertheless, in the present study we do not 
consider any turbulence modeling beyond the scales resolved by the numerical 
method. 

There is much less available information on the snow viscosity and snow rhe- 
ology, in general [DAH + 07]. The viscosity ratio parameter A can be related 
(using equation (23)) to the so-called Schmidt number Sc which represents 
the ratio of the fluid viscosity to mass diffusivity: 



2u 2(1 -5X) 

This number was named after the German engineer Ernst Heinrich Wilhelm 
Schmidt (1892 - 1975) and it is used to characterize fluid flows with simultane- 
ous momentum and mass diffusion processes [YXS02,SS03]. For powder-snow 
avalanches, M. Clement-Rastello reported [CR01] the following values of the 
Schmidt number: 

0.5 < Sc < 1. 



When the Reynolds number is sufficiently high (inertial regime), two scaling 
parameters to be respected are the density ratio 5 = ^+ and the Froude 
number Fr := J° . Recall that in the Boussinesq regime (8 — > 0), only the 

Froude number has to be respected. We quote [PNBNF04] reporting on this 
important issue: 

. . . Satisfying the Froude number and density ratio similarities in the labo- 
ratory means that a very high velocity is necessary, which calls for a very 
large channel. It is not possible to satisfy the density ratio similarity, be- 
cause dimensionless number differs by several orders of magnitude between 
the processes that unfold in nature and those reproduced in the labora- 
tory. . . 

As a result, the interpretation of laboratory results is quite ambiguous. At this 
point, numerical simulation should be considered as a complementary tool to 
physical modeling. 
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3 Numerical methods and simulation results 

In this article we perform Direct Numerical Simulations (DNS) of a snow 
cloud moving down a steep slope. Our solver is based on a freely available 
CFD toolbox OpenFOAM . All computations performed in this study are 3D 
with only one cell in z direction for the sake of efficiency In principle, the 
extension to truly 3D configurations is possible with this code. 

In order to implement model (17) - (19) we modified the standard solver 
twoLiquidMixingFoam which discretizes incompressible two-fluid Navier-Stokes 
equations with Fick's diffusion term in the volume fraction transport equa- 
tion (9). The main modification concerns the momentum balance equation. 
Namely, we had to incorporate two nonconservative terms 2v t 'Vv Vp and 
2uVpVv. For information, we provide here a piece of the code written in 
the internal OpenFOAM language which corresponds to the discretization of 
equation (21): 

f vVectorMatrix UEqn 
( 

fvm: :ddt(rho, U) + 
f vm: : div(rhoPhi , U) - 
f vm: : laplacian(muf , U) - 

U*fvc: : div(DabRho* (mesh. Sf O&f vc : : interpolate (fvc : :grad(rho))))+ 
DabRho*(fvc: :grad(U) O .TO & f vc : : grad(rho) ) - 
DabRho*(fvc: :grad(U) & fvc: :grad(rho)) - 

fvc : : div(muf * (mesh. Sf & fvc : : interpolate (fvc : :grad(U) .T() ) ) ) 

); 

Time derivatives were discretized with the classical implicit Euler scheme. An 
upwind second order finite volume method was employed in space. For more 
details on the retained discretization scheme we refer to [Jas96,Rus02,Ope07]. 
The choice of the finite volume method is justified by its excellent stability and 
local conservation properties (especially in comparison to FEM [ESH04,Eti04,EHS05]). 

3. 1 Description of numerical computations 

In this study we try to reproduce in silico a classical lock-exhange type flow 
with an obstacle placed in the computational domain. The sketch of the fluid 
domain and the initial condition description are given on Figure 4. Experiment 
consists in releasing a heavy fluid which flows under the gravity force along an 
inclined channel. The values of all physical parameters are provided in Table 1, 
while computational domain and discretization parameter are given in Table 
2. 



1(3 




h 



Figure 4. Sketch of the computational domain and initial condition description. 



parameter 


value 


gravity acceleration g, m/s 2 


9.8 


slope, 9 


32° 


friction parameter, a 


0.3 


heavy fluid density, p + , kgjrr? 


20 


light fluid density, p~ , kgjvr? 


1 


heavy fluid kinematic viscosity, v + , m/s 2 


4.8xl0" 4 


light fluid kinematic viscosity, v~ , m/s 2 


l.OxlO" 4 



Table 1 

Values of various parameters used for numerical simulations. The friction parameter 
is used in the partial slip boundary condition (20). 



parameter 


value 


domain height, h, m 


0.8 


total length, £, m 


2.7 


mesh height, Ay, m 


0.0026 


mesh length, Ax, m 


0.0033 


initial avalanche height, ho, m 


0.3 


obstacle height, h s , m 


0.2/i 


obstacle thickness, m 


0.04 



Table 2 

Dimensions of the configuration and mesh. 
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Figure 5. Avalanche at t = 10 s. The color scale ranges from to the maximum 
value 1.0. 

The objective of this study is twofold. Besides presenting numerical results 
on a gravity- driven two-fluid flow, we also would like to shed some light onto 
the interaction process between an avalanche and an obstacle. At the present 
stage we assume the obstacle to be absolutely rigid but this assumption can 
be relaxed in future investigations. Results presented here are obtained for 
the obstacle height h s = 0.2h , where ho is the initial mass height. We would 
like to underline also the presence of the stratification in our initial condition. 
We assume that there is a dense core of the height ho/ 3 composed of the 
pure heavy fluid (0 = 1, consequently p = p + ), which is surrounded by a 
lighter snow layer with the volume fraction <fi = 0.4. All results presented 
below are computed for simplicity with the usual no-slip boundary condition 
v — 0. However, we performed several tests with the partial slip condition (20) 
leading generally to better results with respect to the snow entrainment at the 
bottom. Moreover, the friction coefficient a allows for some tuning depending 
on local soil conditions. 



3.2 Simulation results 

Snapshots of the volume fraction evolution are presented on Figures 5-7. 
At the beginning, the initial rectangular mass gradually transforms under the 
force of gravity into more classical elliptic form (see Figure 5). Then, this 
mass enters into the sliding regime (see Figure 6) until the interaction with 
the obstacle (see Figure 7). 

Our simulations clearly show that a Kelvin-Helmholtz type instability devel- 
ops locally during the propagation stage [Hel68,Kel71,Cha81,DR04]. Previous 
simulations involving Adaptive Mesh Refinement (AMR) techniques confirmed 
this observation [Eti04,EHS05]. On the other hand, the interaction process 
with the obstacle creates a jet directed upward. This jet has a mushroom-like 
shape typical for Rayleigh- Taylor instability [Ray83,Tay50,DR04]. 

Several authors pointed out an intriguing feature of the avalanche type flows 
[DGA01,RH04]. Namely, it was shown by radar measurements that the max- 
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Figure 6. Avalanche at t = 25 s. The color scale ranges from to the maximum 
value 0.984. 



Figure 7. Avalanche at t 
value 0.553. 



60 s. The color scale ranges from to the maximum 




Figure 8. Velocity field magnitude at t = 25 s. The color scale ranges from m/s 
to the maximum value 3.07 m/s. 

imum velocity inside the avalanche exceeds the front velocity by 30% - 40%. 
For this purpose we visualize the velocity field magnitude during the prop- 
agation stage (see Figures 8 and 9). Qualitatively, our computations are in 
conformity with these experimental results. 

The same simulation was also performed with the so-called (in this study) 
standard model used in [ESH04] and implemented also in Open FOAM as the 
twoLiquidMixingFoam solver: 



V-t7=0, 
d t p + v ■ Vp = 2z/Ap, 
pd t v + p(v- V)v + Vp = pg + V- (2pD(v)) 

The standard model differs from governing equations (17) - (19) essentially 
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U Magnitude 
57 




Figure 9. Velocity field magnitude at t 
to the maximum value 6.57 m/s. 



60 s. The color scale ranges from m/s 




Figure 10. Magnitude of two velocity fields difference computed according to the 
standard and new models at t = 60 s. The color scale ranges from m/s to the 
maximum value 2.91 m/s. 




"0.116 

-0.0376 
1-0.191 

-0.344 



Figure 11. Difference of the volume fraction distributions computed according to 
the standard and new models at t = 60 s. The color scale ranges from —0.344 to 
the maximum value 0.269. 

by two non-conservative terms in the momentum balance equation along with 
small differences in the pressure definition. 

The magnitude of the difference between velocity fields obtained with the 
standard and novel models is represented on Figure 10. One can see that 
differences are not negligible and attain its maximum value of 2.91 m/s near 
the lower boundary. Corresponding difference field of two volume fraction 
distributions are shown on Figure 11. Here the maximum value of the difference 
0.344 is attained in the avalanche front. 

During the simulation we also computed the kinetic energy evolution. The 
comparison result is presented on Figure 12. The energy grows almost linearly 
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Figure 12. Kinetic energy evolution during the simulation. The continuous blue line 
corresponds to computations performed with equations (17) - (19), while the dashed 
black line is produced using the standard model. 

in time during the propagation stage for both models. The differences start 
to appear just before the interaction process with the obstacle. These results 
indicate that proposed above model (17) - (19) allows for higher energy levels 
under equal numerical conditions. 

Remark 6 The idea to use the kinetic energy loss to estimate the efficiency 
of a dike was already proposed by Beghin and Closet in 1990 [BC90J. However, 
they had very limited information on the flow structure ( especially velocity and 
density profiles). That is why they decided to approximate this quantity by the 
ratio \U 2 — U' 2 \ /U 2 . Here U' is the front velocity at certain distance below the 
dike and U is the front velocity of the reference avalanche measured at the 
same point. 



3. 3 Impact pressures 



For many practical applications we have to estimate the loading exerted on a 
structure by an avalanche impact. Incidentally, the avalanche hazard level is 
attributed depending on the estimated impact pressure values [Lie06]. More- 
over, this information is crucial for the design of buildings and other structures 
exposed to this natural hazard. 
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In engineering practice, it is common to determine the impact pressures ac- 
cording to the following formula [MS93]: 



where K is a parameter depending on the obstacle configuration, p is the 
average avalanche density and Uf is the front velocity. For small obstacles it is 
advised to take K = 1 and for big ones K = 2 sin a, where a is the incidence 
angle. However, as it is pointed out in [B091], it is difficult to estimate the 
maximum pressure exerted by an avalanche since we have only very limited 
information on the vertical structure of the flow. 

Remark 7 The given above formula (24) is applicable to avalanches in in- 
ertial regime. When we deal with a gravity flow regime (in Fluid Mechanics 
we call it the Stokes flow [HB83]), the situation is more complicated since the 
flow is governed by the rheology which is essentially unknown [AM04J- In this 
case, engineers use another expression [ABB + 06]: 



For aerosol avalanches, Beghin and Closet [BC90] proposed the following em- 
pirical law to estimate the impact pressure: 



where K a (z) is a dimensionless factor taking into account for the velocity 
variations in the upward direction. They also suggested an idealized form of 
the factor K a (z): 



where h is the impacting avalanche height. It was shown later [NB03] that this 
approximation underestimates the dynamic pressure in all parts of the flow. 

The big advantage of the presented here approach is that we have the com- 
plete information on the mass but also the pressure distribution in the whole 
domain. On Figure 13 one can see the distribution of the dynamic pressure 
at t = 60 s during the impact process. It is important to note that an "as- 
piration" zone is revealed near the obstacle base along with a high pressure 
field at the top. In civil engineering this shear-type loading is known to be 
particularly dangerous for structures. 



P d = K Pre{ = KpU 2 f 



(24) 



P d = 2pg(h-z). 





(25) 
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Figure 13. Dynamic pressure distribution in the computational domain at t = 60 s. 
The color scale ranges from —77.5 Pa to the maximum value 4.04 Pa. 

The methodology presented in this study allows to determine the impact pres- 
sures with required accuracy. The pressure profiles along the impacted wall 
can be easily extracted from numerical computations, thus, replacing empir- 
ical formulas such as (25). Our numerical experiments show also that the 
form of the vertical pressure distribution does not change drastically when 
we vary the obstacle height in some reasonable limits. It means that precom- 
puted pressure profiles may be scaled and reused for different obstacles within 
engineering accuracy [BC90,MS93]. 



4 Conclusions and perspectives 

In the present paper we derive a novel model for the flow of two miscible 
inhomogeneous fluids. The mixing between two fluids is assumed to be mod- 
eled by the Fick's law (6). Then, we reformulate the model in terms of the 
so-called volume-based velocity [Bre05a,Bre05b] in which the flow is exactly 
incompressible. It is also shown that the novel model allows for a physically 
sound evolution of the kinetic energy. 

We show some preliminary results on the numerical modeling of powder-snow 
avalanches. We compute the evolution of an avalanche from the beginning 
until hitting against a trapezoidal obstacle. Impact pressures extracted from 
these computations can be used for the design of protecting structures. 

For the moment the obstacle is assumed to be absolutely rigid. In future stud- 
ies this assumption can be relaxed and efforts exerted on the obstacle could be 
coupled with a solid mechanics part [KB01]. This research axis seems to be still 
underexplored. However, there is already interesting experimental material on 
the avalanche interaction with solid obstacles [LSBH95,Pri03,PNBNF04,NB03,BR04]. 

There is also an important question of boundary conditions. All computational 
results presented above were obtained with the usual no-slip condition. The 
development of other types of physically sound boundary conditions which 
lead to well-posed mathematical problems is highly needed. 
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Recall that a powder-snow avalanche front can reach the speed up to 100 m/s. 
Consequently, the Mach number Ma attains relatively high values: 

Ma := ^ « 0.3, 

c s 

where c s is the sound speed in the air. It means that compressible effects 
may be important during the propagation stage. In general, impact events 
are followed by strong compressions. Hence, in future works we are plan- 
ning to take into account the compressibility in some weak and strong senes 
[DDG08b,DDG10,DDG08a,Dut07] and to perform the comparisons with the 
present results. Obviously, at this point a validation against experimental data 
is highly recommended. 

The rheology of avalanches should be further investigated [AM04] and future 
models will take this information into account. At this stage, close collabora- 
tion between physicists and mathematicians is needed to bring the answers on 
challenging questions. 
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